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THE  USE  OP  DIRECT  METHODS  FOR  THE  SOLUTION  OF  THE 
DISCRETE  POISSON  EQUATION  ON  NON-HECTANGULAR  REGIONS* 

J.  Alan  George 

1.  Introduction 

In  recent  years  several  special  direct  methods  have  been  developed 
for  solving  the  discrete  Poisson  equation  on  rectangular  domains.  These 
methods  take  advantage  of  the  regular  block  structure  of  the  coefficient 
matrix,  and  some  of  them  require  an  amount  of  computation  which  is  close 
to  being  directly  proportional  to  the  number  of  grid  points  (equations) 
in  the  discretized  problem.  Dorr  [4]  presents  an  excellent  survey  of 
these  methods.  A  considerable  number  of  these  algorithms  suffer  from 
numerical  instability  and  are  not  suitable  for  large  problems.  An  analysis 
of  stability  of  several  methods  appears  in  [10]. 

In  this  paper  we  describe  ways  in  which  these  direct  methods  can  be 
used  to  solve  non-rectangular  Poisson  problems.  We  will  not  concern  our¬ 
selves  with  which  of  the  direct  methods  is  to  be  utilized;  we  merely  observe 
that  a  number  of  satisfactory  ones  are  available.  Notable  among  them  are 
Buneman 1  s  version  of  the  method  of  odd/even  reduction  [1,2],  and  methods 
based  on  Fourier  analysis  [7,8]* 


This  work  was  supported  in  part  by  the  Office  6f  Naval  Research  under 
grant  N0013-67-A-00112-0029,  the  Atomic  Energy  Conmission  under  grant 
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The  basic  procedure  is  as  follows.  The  domain  R  of  the  given  problem 
is  enclosed  in  a  rectangle  over  which  a  uniform  mesh  is  placed.  The  usual 
five-point  Poisson  difference  operator  is  applied  over  the  entire  rectangle, 
yielding  a  block  tri-diagonal  system  of  equations.  The  given  problem,  how¬ 
ever,  determines  only  those  elements  of  the  right-hand  side  which  lie  in  R; 
the  remaining  elements  can  be  treated  as  parameters.  Furthermore,  the 
"solution"  of  the  enclosing  rectangular  "problem"  which  we  have  generated 
will  have  certain  constraints  imposed  upon  it  by  the  presence  within  the 
rectangle  of  the  boundary  S  of  the  given  (or  imbedded)  problem.  Dirichlet 
boundary  conditions  will  require  the  solution  on  the  rectangle  to  have 
specified  values  at  grid  points  which  lie  on  S;  other  types  of  boundary 
conditions  will  require  specific  relations  to  hold  between  values  at  grid 
points  lying  on  and/or  adjacent  to  S. 

We  now  summarize  our  situation.  We  have  a  fast,  efficient  method  for 
solving  a  specific  system  of  equations,  and  we  cannot  delete  or  modify 
equations  of  the  system  because  the  method  depends  upon  the  structure  of  the 
coefficient  matrix.  We  generate  a  system  of  equations  which  has  this  ap¬ 
propriate  form,  but  for  which  some  of  the  right-hand  sides  are  unspecified,  and 
where  the  solution  must  satisfy  certain  constraints.  This  paper  describes 
methods  for  solving  this  problem. 


2.  Rotation  and,  a  Representative  Problem 


For  definiteness,  we  consider  the  following  problem: 


(2.1) 
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We  superimpose  a  uniform  grid  on  the  rectangle  3,  and  for  simplicity 
we  assume  that  T  lies  on  grid  points  and  on  lines  adjoining  adjacent 
grid  points.  Approximating  the  differential  operator  with  the  usual 
five-point  difference  operator,  and  writing  out  an  equation  for  every 
grid  point  in  S,  we  obtain  an  N  x  N  system  of  equations 

(2.2)  Av  -  h. 


where  the  vectors  v  and  h  are  defined  on  the  grid,  and  N  is  the  total 
number  of  grid  points  in  the  rectangle.  For  expository  purposes  only,  we 
write  (2.2)  in  the  following  partitioned  form; 


(2-3) 
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where  the  vector  partitions  with  subscripts  Q,  R,  and  T  contain  elements 
corresponding  to  grid  points  lying  in  Q,  R,  and  on  T,  respectively. 

We  will  denote  the  number  of  elements  in  these  partitions  by  N^,  NR,  and 
R^.  We  emphasize  that  this  reordering  cannot  be  done  in  practice  because 
the  special  direct  methods  depend  upon  A  having  the  regular  block  structure 
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which  occurs  only  if  the  grid  points  are  numbered  row  by  row  or  column 
by  column. 

It  should  be  clear  that  if  h^  and  h^  are  assigned  values  so  that 
the  solution  to  (2.3)  satisfies  the  boundary  conditions  (i.e.,  if  has 
the  correct  values),  then  v  will  be  the  correct  discrete  solution  to  our 
given  problem.  In  Section  3  a  method  is  presented  for  finding  the  values 
to  be  assigned  to  h^  and  h^  so  that  the  above  is  achieved. 
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3*  Direct  Solution:  Method  1 


3Ms  method  has  been  described  by  Hockney  [  7  ],  and  is  closely  connected 
to  the  discrete  Green's  function  [6]*  Formally,  we  can  invert  the  parti¬ 
tioned  matrix  in  (2.3)  to  obtain 
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,  >  we 

have 

(3-2)  B22  hT  =  vT  -  B21  hj^  -  B23  h^. 

Since  B22  is  non-singular  (it  is  positive  definite),  we  have  set  h^  to 
zero.  We  have  an  efficient  method  for  solving  ( 3 . 1 )  (in  a  reordered  form) 
so  we  Cun  easily  obtain  B21  h^  as  zT  from  the  solution  of  the  system 


(3.3) 
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The  vector  h^,  is  then  obtained  by  solving 

(3**0  Bgg  hT  -  VT  -  *T  • 


Thus,  we  need  B22,  which  means  that  we  need  the  N^,  corresponding  columns 
of  the  inverse  of  the  coefficient  matrix  A.  Ifris  method,  therefore, 
requires  solving  NT  +  2  systems  of  the  form  (2.2),  and  the  solution  of 
the  Nt  linear  equations  (3-^)-  If  we  assume  that  the  number  of  arithmetic 
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operations  required  to  solve  (2.2)  Is  kN,-  then=/  the  total  number  of 

operations  is  about  kfyj,  (N  +  j  N^).  Since  N^,  will  typically  be  Q(/n), 

the  amount  of  work  will  be  roughly  proportional  to  N^.  If  we  suppose  S 

2  U 

is  a  square  with  a  ■  N  ■  10  grid  points  in  it,  and  that  is  2n 
(it  could  easily  be  this  large  for  typical  regions),  then  the  number  of 
arithmetic  operations  required  is  O(n^) .  This  is  not  likely  to  compare 
favourably  with  solving  for  v^  using  SOR,  especially  when  we  consider  how 
little  programming  overhead  there  is  for  the  SCR  process.  Note,  however, 
that  the  matrix  depends  only  on  the  geometry  of  the  problem,  ftius, 
if  we  wish  to  solve  a  time -dependent  problem,  one  with  a  non-linear  right- 
hand  side,  or  many  problems  with  the  same  geometry,  then  this  procedure  may 
very  well  be  the  best  one  to  use.  It  will  almost  certainly  be  the  best  if 
Nt  and  are  small  relative  to  N^. 


Ifce  factor  k  is  actually  a  very  slowly  increasing  function  of  N, 
of  the  form  l  logg  l  a  constant. 
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4.  Direct  Solution;  Method  II 


An  alternate  approach  which  la  acre  general  than  that  of  the  method 
of  section  3  is  the  following:  We  replace  equations  Ag^v^  +  Agg VT  + 

A^v^  *  h^  in  (2.3)  with  the  equations  Ov^  ♦  Iv^  +  Ov^  ■  v^,  by  adding 
a  suitable  correction  to  A.  Obviously,  the  resulting  solution  v  will  have 
the  correct  v^,  regardless  of  the  value  of  h^.  Defining  P  and  G  by 


(4.1) 


(•)■ 


and  denoting  the  coefficient  matrix  of  (2.3)  by  A,  we  can  write  the 
equation 


(4.2) 
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(A  ♦  FG  jv  .  h. 


It  can  be  shown  [  9  ]  that 

(4.4)  (A  +  P0T)_1  -  A’1-  A-1T(I  ♦  GTA_1F)"1QTa";L. 


Thus,  the  procedure  is 

a)  solve  AW  ■  0; 

b)  solve  Ay^  ■  h; 

2)  compute  yg  ■  O^y^  and  Y  ■  I  +  W^T; 
d)  solve  Yy3  -  y2J 
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e)  solve  Ay^  ■  Fy^  j 

f)  compute  v  -  yx  -  y^  . 

Note  that  this  method  is  very  flexible.  It  allows  us  to  replace  any 
equation  by  another  at  the  expense  of  one  solution  of  (2.2). 

The  amount  of  computation  and  storage  required  ia  virtually  the  same 
as  for  method  I.  However,  since  method  II  is  somewhat  more  complicated, 
method  I  seems  preferable  unless  the  increased  generality  provided  by 
method  II  is  necessary. 
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5*  Iterative  Solution  Based  on  Method  I 


We  now  turn  to  potentially  more  efficient  ways  to  utilize  direct 

methods  to  solve  non-rectangular  Poisson  problems.  Our  basic  problem 

is  to  find  a  solution  to  (3A),  and  the  major  expense  in  the  algorithm 

results  from  the  generation  of  Bgg.  Hence,  we  would  like  to  arrive  at 

(k) 

an  iterative  scheme  which  generates  an  (approximate)  solution  h£,  '  with¬ 
out  actually  requiring  Bgg.  First  note  that  for  an  arbitrary  vector 
h£k),  (3.1)  implies 


(5-1) 


(k)  _  ,  (k) 

v'  '  =  z  +  B  h'  ' 
T  T  22  T 


(5*2) 


B  h(k)  -  v(k)-  z  -  w  +  r(k) 
*22  nT  T  ZT  “  W  +  r  ' 


where 


(5*3) 


W  “  VT  "  ZT 


(5-*0 


(k)  (k) 

r  -  vj,  vT. 


(k) 

Here  r'  '  is  the  residual  of  the  linear  system  (3**0  and  is  the  difference 


between  the  solution  on  T  generated  by  h^,  and  the  required  values  vr 
Our  problem  is  equivalent  to  minimizing  the  quadratic  function 


(5*5)  |  h£  B22  hT  -  h£  w, 

where  we  do  not  know  Bg2  but  can  compute  the  gradient  of  cp.  We  are  obviously 

free  to  use  any  of  the  many  iterative  methods  for  solving  a  system  of 

linear  equations  or  minimizing  a  quadratic  function  that  is  bounded  from  below. 


However,  'because  the  residual  (gradient)  calculation  is  expensive 
it  is  natural  to  use  a  relatively  powerful  function  minimizer  or  linear 
equation  solver.  For  example,  we  could  use  the  conjugate  gradient  method, 
or  one  of  the  several  variable  metric  algorithms  which  have  been  developed 
[3,5] .  In  section  8  we  compare  SOR  to  two  of  these  iterative  forms  of 
method  I,  making  use  of  the  conjugate  gradient  method  in  one  and  the 
Davidon-Fletcher- Powell  algorithm  [  5  J  in  the  other.  We  shall  refer  to 
this  class  of  methods  as  Iterative  imbedding  algorithms. 
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Now  the  matrix  In  the  braces  Is 


Hence,  we  have 
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Hie  rate  of  convergence  obviously  depends  critically  on  ||b  q||,  and 

a,  p 

there  appears  to  be  no  easy  way  to  determine  optimal  a  and  p.  For 

a  =  0,  it  is  easy  to  show  that  p  should  be  set  to  2/(^max+  ^min)  where 

X  and  X  .  are  the  largest  and  smallest  eigenvalues  of  B_~.  B„  _  then 
max  min  0  °  22  0,p 

has  as  its  largest  (in  magnitude)  eigenvalue  |i  ■  (X  -  X  .  )/(x  +  X  .  ), 

max  mm  max  mm 

and  the  iteration  then  converges.  Hie  problem  is,  of  course,  the  difficulty 

in  determining  estimates  of  X  .  and  X  .  .  Numerical  experiments  and 

max  min 

further  analysis  of  this  method  are  currently  being  pursued. 
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One  of  the  most  difficult  problems  in  the  application  of  an  iterative 


process  is  the  determination  of  a  safe  and  meaningful  convergence  criterion. 
For  a  short  and  very  good  account  of  this  problem  with  SQR  see  [8  ]. 

Briefly,  the  problem  is  as  follows:  Since  we  do  not  know  the  true  (discrete) 
solution,  the  error  at  each  stage  of  the  iteration  must  be  estimated  on  the 
basis  of  such  measurable  quantities  as  the  size  of  the  residuals  or  the  size 
of  the  last  correction  vector.  Unfortunately,  small  residuals  or  small 
changes  in  successive  Iterates  do  not  guarantee  correspondingly  small  errors 
in  the  computed  solution.  For  rather  ordinary  problems  the  error  can  be 
several  orders  of  magnitude  larger. 

The  iterative  imbedding  algorithms  seem  particularly  attractive  with 
regard  to  the  above  problem,  as  the  following  theorem  demonstrates. 

Theorem  1.  Let  v  be  the  true  discrete  solution  on  the  enclosing 
rectangle,  and  let  v*  be  the  computed  solution,  where  v*  satisfies 
the  (DLrichlet)  boundary  condition  to  within  some  value  e,  i.e., 


(7-1)  l|vT-  v*|L  <  e. 

Then 

(7.2)  .  ||vR-  vg||  <  €  • 

Proof:  Let  L  be  the  discrete  Laplacian  operator.  Then  the  following 
equations  are  satisfied: 


(7-4) 


Lvr  -  "a, 

1vr-  V 
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Denoting  the  error  in  the  computed  solution  by  e,  we  have  from  ( 7-1) , 
(7.3),  and  (7.4)  that 

(7-5)  LeR  -  ° 

and 

(7*6)  ||eT||w<€. 

Since  -L  is  an  operator  of  "positive  type,"  we  can  apply  the  well-known 
maximum  principle  to  conclude  from  (7*5)  and  (7-6)  that  lie^liao“!!v^*vpll  S  e' 
Thus,  we  can  determine  when  to  stop  the  iteration  simply  by  examining 
the  largest  element  of  je^l  .  Since  it  is  difficult  to  imagine  an  iterative 
scheme  which  would  not  make  use  of  e^  (it  is  the  residual  of  (3*4)),  the 
cost  of  determining  HeJI^  should  be  negligible. 
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8.  Numerical  Experiments 

We  now  present  some  numerical  experiments  for  a  problem  of  the  type 
(2.1),  where  S  is  covered  with  a  square  mesh  having  m  rows,  n  columns, 
and  mesh  width  h,  and  where  Q  is  a  1th  x  ih  rectangle.  Ibe  "southwest" 
corners  of  S  and  Q  are  at  grid  positions  (0,0)  and  ( )  respectively. 

The  implementation  of  the  SQR  algorithm  provides  for  an  initial  ap¬ 
proximate  solution  on  a  coarse  grid  (with  mesh  width  2h)  which  is  then  used 
to  furnish  a  starting  solution  on  the  fine  net  by  using  linear  interpola¬ 
tion.  Biirty  iterations  were  carried  out  on  the  coarse  mesh  to  obtain  the 
initial  solution,  and  these  iterations  and  the  time  required  for  them  are 
not  included  in  the  tables  below.  An  acceleration  parameter  uj  of  1.8 
was  used  on  the  coarse  mesh  for  the  first  25  iterations,  followed  by  5 
iterations  with  w  »  1  to  estimate  the  optimal  uj  ■  w*  for  the  coarse  mesh. 
The  value  tu*  +  -55  (2-iu*)  was  found  to  be  near  optimal. for  the  fine  mesh. 

The  number  of  iterations  reported  for  the  iterative  imbedding  methods 
requires  some : discussion.  Obviously  each  iteration  requires  substantially 
more  work  than  an  SOR  iteration.  The  ratio  will  depend  on  the  size  of  the 
mesh  since  the  computation  required  for  the  direct  methods  is  not  quite 
directly  proportional  to  mn.  Also,  the  relative  sizes  of  N_  and 

a 

N_  +  N  +  N_  will  affect  the  ratio  because  the  SOR  iterations  will  (at 
A  T  ^ 

least  ideally)  only  involve  grid  points  in  R.  A  factor  of  about  10  seems 
reasonable  for  typical  problems  having  fewer  than  20,000  points. 

Ibe  time  required  to  compute  the  right-hand  sides  of  the  equations 
is  not  included  in  the  tables.  All  times  are  for  execution  of  ALGOL  W 
programs  on  an  IBM  360/ 67  • 
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Case  I:  f  ■  (2-100  y^)  cos  (lOx) 
o 

u  ■  g  ■  y  cos  (lOx) 
h  -  0.0125,  m  ■  49,  n  -  127 
(J^)  -  (64,32),  k  -  i  -  10 

Case  II:  Same  as  Case  I  except  (j^Jg)  ■  (20, 40)  and  k  ■  l  ■  20. 


Method 

Iterations 

Maximum 

Error 

Time 

(Seconds) 

S0R 

70 

4. 7x10“ 4 

42 

Imbedding  I 

5 

2xl0“4 

24.0 

Case  I  T  -w  xt 

Imbedding  II 

5 

2x10" 4 

23.7 

Direct 

N.A. 

2x10" 11 

* 

9-6 

S0R 

70 

4.2x10”^ 

41 

Imbedding  I 

6 

2x10" ^ 

28.6 

Case  II 

Imbedding  II 

6 

2xl0‘4 

28.5 

Direct 

N.A. 

2xl0“4 

9.6* 

Imbedding  I  -  method  of  Section  5  using  the  DavicLon-Fletcher-Powell 
algorithm  [  5 ] . 

Imbedding  II  -  method  of  Section  5  using  the  conjugate  gradient  algorithm. 
Direct  -  method  of  Section  3* 

The  maximum  errors  for  the  direct  method  and  the  imbedding  methods 
are  all  the  same  because  the  error  is  due  entirely  to  the  truncation  error 
of  the  difference  operator.  The  error  in  the  discrete  solution  for  these 
methods  is  below  that  level. 


Does  not  include  the  time  required  (approximately  3  minutes  and  6 
minutes,  respectively,  for  Cases  I  and  II )  to  generate  and  decompose 

B22  (see  Seotion  3)* 


16 


9.  Remarks  and  Conclusions 

The  reported  times  at  first  do  not  appear  particularly  impressive, 
although  the  times  required  for  the  imbedding  methods  are  substantially 
less  than  for  the  SOR  process.  It  is  important  to  keep  in  mind,  however, 
that  during  the  calculation  using  the  methods  of  Section  5>  we  have  precise 
information  concerning  how  close  our  computed  solution  is  to  the  true 
discrete  solution.  This  is  obviously  highly  important  in  a  practical 
situation  where  the  solution  to  our  problem  is  not  known.  As  we  mentioned 
in  Section  7,  it  is  extremely  difficult  when  applying  SOR  to  ascertain 
how  close  the  computed  solution  is  to  the  true  discrete  solution.  (For 
example,  the  maximum  change  for  the  last  step  of  SOR  in  Case  I  above  was 
8.1xlO~5.) 

As  one  might  expect,  the  rate  of  convergence  of  the  iterative  imbed¬ 
ding  algorithms  depends  on  N^,.  However,  quite  extensive  experiments  seem 
to  suggest  that  the  number  of  iterations  does  not  increase  very  rapidly  with 
N^,,  and  iterations  are  usually  sufficient.-  Problems  with  singularities 

also  do  not  appear  to  greatly  affect  the  rate  of  convergence. 

When  N_  and  N_  are  relatively  large,  and  R  can  be  subdivided 
W  i 

into  a  number  of  rectangular  blocks  (R  might  be  H-  or  L-shaped,  for 
example),  a  direct  method  described  in  [  2  ]  may  be  more  efficient  than  the 
method  described  in  Section  3.  It  is  not  obvious  if  or  when  its  iterative 
analog  converges  and,  even  if  it  does,  no  a  posteriori  bounds  are  available 
because  the  "parameters"  are  grid  values  lying  in  R  rather  than  on  a 
boundary. 


We  have  not  discussed  the  direct  method  used  to  solve  our  rectang¬ 
ular  problems.  As  ve  mentioned  earlier,  many  of  the  methods  discussed 
by  Dorr  [k]  suffer  from  numerical  instability  and  are  not  suitable  for 
large  problems.  We  have  used  a  method  due  to  Buneman  [1]  which  appears 
to  be  stable  even  for  very  large  problems.  For  a  qualitative  discussion 
explaining  this  stability,  see  [2].  Hockney's  algorithm  POT  1  [7]  could 
in  theory  reduce  the  times  for  the  imbedding  algorithms  and  the  direct 
method  by  a  factor  of  two,  although  in  practice  program  overhead  would 
reduce  some  of  the  advantage  of  the  lower  operation  count. 

.  Note  that  no  use  has  been  made  of  the  particular  geometry  of  the 
problem  we  have  discussed  other  than  it  is  enclosed  by  a  rectangle.  The 
methods  we  have  described  are  applicable  to  arbitrary  domains,  and  their 
efficiency  will  depend  upon  the  subjective  factors  discussed  at  the  end 
of  Section  3. 
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